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Abstract 

We extend the work of Tanase-Nicola and Kurchan on the structure of diffusion processes and 
the associated supersymmetry algebra by examining the responses of a simple statistical system 
to external disturbances of various kinds. We consider both the stochastic differential equations 
(SDEs) for the process and the associated diffusion equation. The influence of the disturbances 
can be understood by augmenting the original SDE with an equation for slave variables. The 
evolution of the slave variables describes the behaviour of line elements carried along in the 
stochastic flow. These line elements together with the associated surface and volume elements 
constructed from them provide the basis of the supersymmetry properties of the theory. For ease 
of visualisation, and in order to emphasise a helpful electromagnetic analogy we work in three 
dimensions. The results are all generalisable to higher dimensions and can be specialised to 
one and two dimensions. The electromagnetic analogy is a useful starting point for calculating 
asymptotic results at low temperature that can be compared with direct numerical evaluations. 
We also examine the problems that arise in a direct numerical simulation of the stochastic 
equation together with the slave equations. We pay special attention to the dependence of the 
slave variable statistics on temperature. We identify in specific models the critical temperature 
below which the slave variable distribution ceases to have a variance and consider the effect on 
estimates of susceptibilities. 
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1 Introduction 



Stochastic processes have applications in many areas of physics. They are a natural way of 
describing the effect of thermal interactions on dynamical systems subject to force fields of 
various kinds. These force fields could be, for example, the inter-atomic potentials represented 
by an energy landscape appropriate to a molecular model [1]. In another type of application, 
stochastic processes provide a method for achieving a desired probability distribution for the 
system. This method is the basis for the theory of stochastic quantisation in quantum field 
theory. It has been exploited in the numerical simulation of quantum field theories [2, 3, 4]. 

One way of investigating the structure of a dynamical system is to subject it to external 
disturbances. This can be accommodated in the stochastic dynamics by inducing small changes 
in the motion of the particle either by changing its initial position or by altering the force field 
guiding the particle throughout its motion. If these external influences are infinitesimal the 
result, for a given sample from the noise ensemble, is an infinitesimal shift in the path of the 
particle. It can be regarded as an infinitesimal line element carried by the particle along the 
path of the original motion. Because the equation describing the evolution of this line element 
is derived from the original stochastic equation and turns out not to contain the noise term we 
refer to it as a slave equation and to the infinitesimal line element as a slave variable. 

The response of the mean position of the system to external influences, that is its suscepti- 
bility, can be calculated once the statistical properties of the slave variables are known. Some 
aspects of slave variable statistics in a one dimensional case have been investigated in [5] . In 
higher dimensions there exists also the possibility of examining the evolution and statistical 
properties of infinitesimal area and volume elements that can be constructed from the the line 
elements. The evolution of area and volume elements plays a significant role in turbulent fluid 
flow and in magnetohydrodynamics [6, 7, 8, 9, 10, 11, 12]. In the models studied in this paper 
these higher dimensional area and volume elements can be thought of as the images of area and 
volume elements in the parameter space of the external fields. Their statistical properties thus 
give rise to higher dimensional susceptibilities. These generalised susceptibilities can be readily 
computed in numerical simulations and their evaluation can help to elucidate the topography of 
the landscape function. In addition we will find that the hierarchy of generalised susceptibilities 
fits naturally into the supersymmetric structure underlying the stochastic dynamical system. In 
turn the supersymmetric analysis makes possible the calculation of these quantities. 

A comprehensive review of material that is relevant particularly to one-dimensional systems 
and the supersymmetry approach we adopt here is contained in reference [13]. Wittcn [14] 
has explained the relationship between supersymmetric quantum mechanics and the topology 
of manifolds. Making use of the idea that diffusion processes can be regarded as quantum 
mechanics in imaginary time, Tanase-Nicola and Kurchan [15, 16] have developed these ideas 
emphasising the diffusion point of view and using it to construct practical simulation techniques 
for picking out the saddle point structure of the landscape function. 

In this paper we extend the work of Tanase-Nicola and Kurchan [15, 16] and show it can 
be developed to compute significant physical observables. We do this in part by developing 
an electromagnetic analogy based in the Principle of Minimum Dissipation. This allows us 
to understand the structure of the physical states of the system in a way that is particularly 
revealing in the low temperature limit. It can also be the basis of an optimised approximation 
scheme at higher temperatures though we do not pursue this idea here. 

The structure of the paper is as follows. In section 2 we describe the basic stochastic model 
that we study. In section 3 we derive the slave equations, first for an infinitesimal change 
in the initial position of the particle and second for a weak external field. We show how they 
can be encoded in the supersymmetric formalism by introducing an appropriate supersymmetric 
Hamiltonian. Section 4 contains a brief exposition of the numerical method we use for simulating 
the stochastic differential and slave equations. In section 5 we give, for completeness and to 
establish notation, a brief account of the state space on which the Hamiltonian acts and an 
outline of the structure of its eigenstates. In section 6, we develop an electromagnetic analogy 
that helps us to understand the structure of the eigenstates of the supersymmetric Hamiltonian 
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and their time evolution under various circumstances. The analogy is developed in section 7 
where we show that the evolution of the eigenstates obeys a Principle of Minimum Dissipation. 

In section 8, the low temperature limit in which it is possible to compute low lying eigen- 
values of the supersymmetric Hamiltonian and associated decay exponents is discussed. The 
application of the supersymmetric theory to the case of an external electric field together with 
a derivation of the associated Einstein Relation is set out in section 9. Explicit examples of the 
calculations, in one and two dimensions, are exhibited in sections 10 and 11, where we compare 
the results of various numerical techniques. Finally we we discuss the results and their potential 
significance in section 12. 



2 Basic Stochastic Model 

In this paper we consider a particle with three degrees of freedom, {A} — {Ai, A 2 , A 3 }, that 
obeys the stochastic differential equation of the form 

X n = -<j> n (X) + W n (t) , (1) 

where W n (t) is a stochastic process for which 

(W n (t)W m (t')) = 5 nm G(t - t') , (2) 

for some correlator G(t — t') . The angle brackets denote averaging over the noise ensemble. 
Aspects of our discussion are independent of the precise nature of the noise W n (t) . Ultimately 
however we will be interested in the case of white noise. From now on therefore we will assume 
that 

(W n (t)W m (t')) = 2T5 nm S(t-t') . (3) 

We will refer to T as the temperature of the system. It is implicit in our assumptions that the 
ensemble of white noise processes is rotationally invariant in the particle configuration space. 
While this is not the most general case it is sufficient for illustrating our approach. 

The drift term in eq(l) is assumed, as indicated above, to be derived from an energy land- 
scape function <j>(X) . That is (f> n (X) = d n <fr(X) . Higher derivatives of <fi(X) will also be 
indicated with suffixes, thus 4> mn (X) = d m d n (/)(X), and so on. It is convenient and not really 
restrictive to assume that </>(A) rises like a power of \X\ as \X\ — » oo . As a consequence the 
particle will remain confined to an essentially finite region in A-spacc throughout its motion. If 
we assume that at time t = the particle is at the point {y n } then at later time t, for a given 
element {r] n (t)} of the noise ensemble, the motion of the particle will be that solution, {X n (t)}, 
of eq(l) that satisfies {A„(0)} = {y n } ■ The probability distribution for the position is 

P(x,t) = (6(x-X(t))) . (4) 

It is standard theory that P(x, t) satisfies 

^P{x,t)=Td n {d n + pif> n {x))P{x,t) , (5) 

where (3 = l/T . Our assumption about the form of (j)(X) implies that as t — > oo P(x, t) — > Po(x) 
where 

Td n (d n + 0(f> n (x))Po(x) = . (6) 
The static solution is unique and has the form 

P {x) = ^e-MW , (7) 

where 

Z = [ dxe~ mx) . (8) 
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Here for compactness we use dx = dx\dx<idx% . Clearly 



(d n +f3<f> n (x))P (x)=0 . (9) 

It is also convenient to introduce 

PoW^i^j^e-^W . (10) 

We have 

(d n + \<t>n{x)) P0 {x) = . (11) 

Eq(5) incorporates the operator H where 

H = -d n (d n + p(P n (x)) ■ (12) 
It is convenient to construct an associated hermitian operator H where 

Ho = (poWr'HoPoix) = -(d n - l -(3cj> n {x)){d n + l -(3<l> n {x)) . (13) 

The two operators H and H have the same eigenvalues but the latter has the additional 
property that its left and right eigenfunctions are identical. We note that 

H = -d 2 n + ^(3 2 ( ( f >n (x)) 2 -^(3<j )nn , (14) 

and finally that eq(ll) implies that 

H oPa (x) = . (15) 

3 Slave Equations 

For a number of reasons both mathematical and physical it interesting to examine the effect 
on the motion of the particle of small external changes. In this paper we consider two types of 
change, namely, a change of initial conditions where {y n } {yn + ^Un}, and the application 
of an external field where the drift velocity is modified by making the replacement <j) n (X) — > 
cf) n {X) — J a fna{X). The functions {f na {X)} are a basis set for the external influences in which 
we are interested. A significant case is one where f na {X) = d n f a (X). This is equivalent to the 
replacement^^) — > (j>(X) — J a f a (X). The choice f a (X) = X a , a = 1, 2, 3 describes the action 
of an electric field on the motion of the particle. We will be interested also in situations in which 
fna(X) is not the gradient of a scalar field. An external magnetic field falls into this category. 

We assume that all the effects are weak so it is sufficient to consider a limit in which for a 
given noise function the induced change, X n — > X n + 5X n , is such that SX n = O(Sy) or O(J). 

3.1 Variation of Initial Position 

The effect of an infinitesimal variation of the initial position of the particle can be computed 
from the behaviour of the quantities {S m „} where 

" <9A m 

^mn = • (16) 

dy n 

The differentiation is performed for a fixed sample {W„(t)} of the noise. These variables satisfy 
the equation 

^mn = — 4>mk(X)^kn ■ (17) 

We refer to the {S m „} as slave variables and to eq(17) as a slave equation. The slaves depend on 
the noise only through the original position variables {A m } but the latter are not influenced by 
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the evolution of the slaves. These variables were introduced by Graham [17] who first elucidated 
their relationship with the underlying supersymmetric structure of the theory. He identified the 
relationship of the associated decay exponent with a supersymmetric eigenvalue problem. 

Because it is a useful introduction to the properties of slave variables we will briefly review in 
our notation some of his results. The first step is to introduce the joint probability distribution 
P(x,£,t) where 

P(x,t,t) = (8(x-X(t))5(£-E(t))) . (18) 
It is standard theory, given eq(l) and eq(17), that 

d d 

7 -P(x,Z,t)=Td n (d n + P<t> n )P(x,Z,t) + 7 —(t>rnk(x)ZknP(x,U) • (19) 

We introduce the quantities P mn (x), where 

P mn (x,t) = Jd^ mn P(x,U) • (20) 

At time t 

{^mn) = ^mn J dxd£^ mn P t) / dxP rnn ( K X 7 t^) . (21) 

Here is the volume element in £-space. From eq(19) we can show that 
d 

—P mn (x,t)=Td n (d n + P(f) n )P mn (x,t)-(t) mk (x)P kn {x 1 t) . (22) 

Following Graham [17], we introduce a set, {a n }, of anticommuting variables, one for each 
of the variables {x n }, together with their hermitian conjugates, {a^}, and a fermionic ground 
state |0) . The new fermionic variables satisfy the (Clifford) algebra 

{a m ,a n } = {a^al} = , (23) 

and 

{a m ,<4} = 5 mn . (24) 

The fermionic ground state satisfies 

a m \0) =0 , (25) 

and is normalised so that (0|0) = 1 . 

We now augment the operator Hq by defining H where 

H = H + P(j) m n{x)a) m a n , (26) 

and define the state ip so that 

Vo=Po(z)|0) • (27) 

It is easily checked that 

Hil> = . (28) 

We now introduce F mn (x,t) where 

Pmn{x,t) = F mn (x,t)P (x) , (29) 

and define the state 

tpn(t) = F mn (a;,i)a] n Vo • (30) 
Eq(22) can now be put in the form 

j^„(t) = -THi> n (t) . (31) 
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The formal solution to this equation is 

il> n {t) = e- TH V™(0) . (32) 

It follows that the long time behaviour of F mn (x,t) and hence that of P mn (x,t) and finally of 
£„ m is dominated by the lowest eigenvalue, E, of H in the one fermion sector of the state space. 
Typically we expect that at large times 

£,mn ~ 6 ■ (33) 

Therefore the decay exponent v, that governs the large time statistical evolution of line elements 
is given by 

v = TE . (34) 

This is the essence of Graham's results [17] . We will examine later the supersymmetric structure 
of H. 

The above approach can be extended to area elements. We consider the infinitesimal area 
element at t = which has the form dy m A dy n . At a later time t it has evolved into the element 



That is 
where 

Clearly 
where 



dX m A dX n = E mk ^ n idy k A dyi . (35) 



dX m A dX n = E mnM dy k A dyi , (36) 



(38) 



^mn.kl 2 (^<7rik£,nl Cm/Cnfc) ■ (39) 



We define P mn ,ki{x,t) and F mnM (x,t) so that 

Pmn.kl{x,t) = F mnM {x,t)P {x) = j ' d€€ mn ,klP(x,€,t) . (40) 

From cq(19) we can show that 
d 

g^Pmn,kl( X ^) = T dn{d n + (3(f) n {x))P mnM (x,t) 

- (t) m k'(x)Pk'n,kl( x ^) - <i>nk'{x)P m k',kl{x,t) . (41) 

We define 

Ipkl = Fmn,kl( x > t ) a ln a li'&0 ■ (42) 

Eq(41) is equivalent to the result 

^ki = ~TH^ kl . (43) 

It follows that the behaviour of ip k i at large times, and hence that of P m n,ki(x,t) and (E m n,kl)> 
is dominated by the lowest eigenvalue of H in the two-fermion sector. This is again one of 
Graham's results [17]. 

Clearly the argument can be generalised to higher dimensional volume elements. We will 
not exhibit the argument in detail. The important point is that the time evolution of a state, 
tp, at any fermionic level is governed by the equation 

^ = -THi> . (44) 

Decay exponents for differential volume elements of any dimension can be calculated through 
the eigenvalues of the supersymmetric Hamiltonian. 
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3.2 External Fields 



The effect of an external field can be induced by replacing eq(l) with 

X„ = —4> n (X) + J a f na (X) + W n (t) . (45) 

Of particular interest is the electric case where f na (X) — d n f a (X). We can think of the set 
{f a (X)} as a set of electric potential functions. If we make the further specialisation 

f a (X)=X a , (a -1,2, 3) , (46) 

the resulting field is an uniform electric field (J 1; J 2 , J3) . The magnetic case is one in which 
d m fna(X) — d n fma(X) ^ and fna{X) is not the derivative of a function {f a (X)}. We will 
consider the magnetic case in detail later. The following derivation covers both cases. 

If we solve eq(45) for a particular sample of the noise process with a given initial condition, 
X(0) = 0, say, then at a later time, t, the value of X will depend on the value of J„ and we 
can regard the solution at this time as providing a map from J-space to X-space. Therefore 
curves and surfaces in J-space will have images in X-space that are also curves and surfaces. 
We can examine the differential version of this map by introducing appropriate slave variables 
E mn given by 

~ma • [*<) 

we have then 

dX m = E ma dJ a . (48) 

By differentiation cq(45) with respect to J a we can obtain an equation for the evolution of S 
namely 

S ma = -(f> m k(X)E ka + f m a(X) , (49) 

Of course we have set J a = after differentiation. The initial condition for the slave is S ma (0) = 
since we assume that X n (0) is independent of J a . 
Eq(48) implies that 

d(X m ) = (dX m ) = (E ma )dJ a . (50) 

that is 

JL(X m ) = (~ ma ) . (51) 

Thus (S m o) can be understood as a set of susceptibilities of the system at time t with respect 
to the various contributions to the external field in the weak field limit J a — > . 

Eq(l) and eq(49) imply that the joint probability distribution, P(x,£, t), given by 

P(x,U) = (S(x-X(t))8(t-Z(t))) (52) 

satisfies 

-P{x, t) = Td n (d n + M n (x))P{x, & t) 

d 

+ ^7 {4>mk£,ka{x) - f ma (x))P{x,£,t) . (53) 

An important distinction between this case and the previous one is that we expect the joint 
probability distribution to approach a static form at large time. That is P(x,£, t) — > P(x,£) as 
t — > oo, where 

d„(d n +0<f> n (x))P(x,O +l3 7 ^—^ mk (x)^a - fma(x))P(x,Z) = . (54) 
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We shall be mainly interested in this asymptotic static situation from now on. Following the 
pattern established in subsection 3.1 we define, for this new case, P ma (x) and F ma (x) so that 

Pma(x) = F ma { X )P {x) = j dj& m aP(x,Q . (55) 

Hence 

(H mo ) = £ ma = J dxF ma (x)P (x) . (56) 

We can interpret F ma (x) as the contribution to the susceptibility from the particles in the 
neighbourhood of the point x. From eq(54) we can deduce that 

d k (dk + PMx))Pma{x) - f3<f>mk{x)P k a(x) + (3 f ma {x)P (x) = . (57) 

Invoking the fermionic variables introduced above and setting 

"0a = F ma (x)al n tpo , (58) 
we find that eq(57) can be expressed in the form 

H$ a -Pf ma {x)a\j> = Q . (59) 

3.2.1 Area Elements 

The discussion can be extended to area elements. At time t the solutions of eq(45) map the 
J-space area element d.J a A dJb into the X-space area element dX m A dX n thus 

dX m A dX n = E ma E n bdJ a A dJb — Z-mn,abdJ a A dJb i (60) 

where 

"mn,o6 ^{^ma^nb ^na^mb) • (61) 

We can think of (S mni0 &) at large t as a generalised area susceptibility. The relationship of this 
susceptibility for the supersymmetric structure can be understood by defining P mn ,ab{x) and 
F mn ,ab(x) so that 



Pmn,ab{x) — F mna b{x)P^{x) J d^ rnna bP{x i ^) , 



(62) 



where 

£,mn,ab — ^{Cma^nb £,na^mb) • (63) 

Hence 

(S mn , a6 ) = imn,ab = j dxF mn , ab {x)P {x) . (64) 

From eq(54) we can show that 

d P {d p + P4> p {x))P mn , ab (x) 

fi^rnm' -^m'n,a&(^) $$nn' ' (^)-^mn' ,a&(^) 
+ P(fma(x)P n b(x) + fnb{x)P ma {x) - f m b{x)P n a(x) - fna{x)P mb (x)) = . (65) 

If we now define 

"0a6 = Frnn.ab^O^a^Q , (66) 

then eq(65) can be put in the form 

Hip ab - ^(f ma (x)al n i() b - f nb {x)alip a ) = . (67) 

From eq(67) we can calculate the area susceptibility in terms of the eigenstates of H at the 
second fermionic level. The procedure can be generalised to deal with volume elements (of 
any dimension). Higher order susceptibilities of this kind fit naturally into the supersymmetric 
approach to the diffusion model. However, while we feel that they are an interesting area 
for study, in contrast to the standard susceptibility to the action of an external field, their 
observability is unclear. 
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4 Simulation 



We will test some of the theoretical predictions of the stochastic model by simulation. The first 
order algorithm corresponding to eq(l) with correlator given by eq(3) is a process in which for 
a time step 5t, X — > X + 5X where 

SX n = -<t> n (X) St + VZTSi Ti n . (68) 

Here r] n is a Gaussian random variable with correlator 

(VmVn) = $mn ■ (69) 

The errors in the simulated probability distributions are 0(6t) . In order to control errors more 
easily we use a second order algorithm with an intermediate step X — » X^ = X + SX^ where 

SXW=-^4> n (X)6t + Vf5tTi^ . (70) 

The final step is 

5X n = -<j )n (X^)St + VWt(r i P+r l W) ' (71) 

where r]n^ and rj n 2 ^ are Gaussian random variables of the same type as r\ above. This procedure, 
which is a stochastic second order Runge-Kutta technique [] can be extended to include the 
slave equation in an obvious way. The resulting probability distributions have errors that are 
0(St 2 ) . 



5 Supersymmetric Formulation 

Apart from the straightforward case in which <p{X) is a quadratic function, the equations for X 
and S are intrinsically non-linear. It is a remarkable and significant fact that associated with such 
non-linear processes is an underlying linear algebraic structure [14, 15, 16, 18, 17, 19, 20]. The 
supersymmetric formulation of the theory is a powerful way of organising this linear algebraic 
structure. Although well known, see in particular [15, 16], in the interests of completeness and 
to establish some notation we give an account of it here. 

The supersymmetric structure is made evident by the construction of the operators 

Q = -iai (di - \p<t>i\ , (72) 



2" 



and 



Q f = -iaj + ^8&) • (73) 
They satisfy 

Q 2 = (Q t ) 2 = . (74) 
The combinations Q\ = Q + and Q 2 = i(Q — Q^) are hermitian and satisfy 

Q1Q2 + Q2Q1 =0 . (75) 

We also have 

Q\ = Ql = H , (76) 

where 

H = QQ^+Q^Q , (77) 

and H is given by eq(26). 
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5.1 The Kernel of H 



From eq(26) we see that ipo is in the kernel of H . Eq(76) then tells us that ^ is in the kernel 
of Q\ and that of Qi also. In fact any eigenstate in the kernel of H is annihilated by both Q\ 
and Q2 and hence by both Q and . 

The operator A = a\a,k counts the number of a-excitations in a state. It commutes with H 
and can be used to classify the eigenstates of H. In particular the basis states in the kernel of 
H can be chosen so that they are each associated with an eigenvalue of A, that in our three- 
dimensional model can take one of the values 0, 1, 2, or 3. If there is only one basis state in the 
kernel of H then it can be taken to be Vo and 

A^ o = . (78) 

In fact because our model involves only a simple flat euclidean space, the kernel of H is one 
dimensional and spanned by t}> [14]. 

5.2 Eigensubspace structure of H 

We also see that if \ is an eigenstate of H with eigenvalue E ^ then so is the orthogonal state 
x' — QiX ( or equivalently QiX)- That is, the eigenstates of H that are not in its kernel come 
in pairs with equal eigenvalues. There is no corresponding constraint on the eigenstates lying 
within the kernel. 

On any subspace W\ with A = A and < A < 3, that does not intersect the kernel of H we 
have the projection operators 

0(-)=Qt Q A : , (79) 

and 

tt (+) =QQ f 4 • (80) 

They satisfy 

fi(+) = 1 , (81) 

and other appropriate relations. It follows that, in general, W\ is a direct sum of two orthogonal 
subspaces = Q^Hx that satisfy 

n^n[ ±] = H[ ±] , and Q^H x t) = . (82) 

Therefore the basis states of Ti.\ can be chosen so that they lie in either H x + ^ or % x ^ . In an 
obvious notation we have 

H X Xa - A Aa X\ a ■ [od) 
It is then obvious from the discussion that we could choose phases so that 

x[-J = r J r _ Qt x W , (84) 



\f E (X-l)a 



with E{J = E$l 1)a and 



& nd E ( £ = E ( - X l 1)a . 



= -jL^Qx^ , (85) 



J Xa 



It follows then that any state in %\ of the form Q^x lies in • ^ ^ s a ^ so * rue tnat an y 

state that lies in H x ' can be put in the form where \ lies i n T~L(\-i) ■ The state \ is 
not unique but the ambiguity can be removed by requiring that \ li es i n ^(A-i) ■ That is, by 
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imposing the condition Qx = . Similarly a state that lies in H x can be put in the form Qp 
for some state p that lies in and that satisfies p = . 

The above analysis has to be modified for the extreme subspaces Ho and V.3 . The subspace 
Ho is a direct sum of H { +) and the one dimensional subspace containing ip ■ The subspace H.3 
is identical with . 



6 Electromagnetic Analogy 

An investigation of the eigenstates and eigenvalues of the supersymmetric Hamiltonian is equiv- 
alent to studying the evolution equation, eq(44) at each fermionic level. It is helpful, then, to 
think in terms of an analogy with electromagnetic systems. The analogy makes it natural to 
introduce the Principle of Minimum Dissipation that we establish in section 7. In turn this will 
help us to establish in a natural way the structure of the eigenstates of the supersymmetric 
Hamiltonian and in particular, in the low temperature limit as elucidated by Tanase-Nicola and 
Kurchan [15] 



6.1 Zero Fermion Level 

Consider, for example, a state ip in the zero fermion sector. It has the form 

1> = V(x)ifo ■ (86) 

In the context of our proposed electromagnetic analogy, it is convenient to associate with it a 
charge density distribution p(x) where 

p{x) = V{x)P {x) . (87) 

(Of course p and V may also have a dependence on t that we leave implicit). 

The image of ip in the one fermion sector is x = ■ It follows immediately that Q^x = 

and therefore that x^T~L\ ' • Evaluating x explicitly we find 

X = ai{d n + 1 -^ n {x))V{x)i,o = {d n V{x))aii,o ■ (88) 

We continue the electromagnetic analogy by associating with x an electric current density j 
where 

i = -TPo(x)VV(x) . (89) 

The justification for these identifications is that, if we assume that satisfies eq(44), we find 
that 

|p=-V-j , (90) 

which is indeed the law of conservation of charge. We can complete the analogy by identifying 
V(x) as the voltage distribution and —\7V(x) as the electric field acting on the charge. The 
static probability density distribution Po(x) plays a dual role. It is the capacity per unit volume 
and also determines the local conductivity, which is TP (x) . 

The eigenvalue equation for the supersymmetric Hamiltonian in the zero fermion sector can 
be re-expressed in the analogy by making the replacements V = e~ TEt , p = p^e~ TEt and 
j — > j^e~ TEt , where the superfix (0) indicates that the quantity is evaluated at t = and E is 
the eigenvalue of the supersymmetric Hamiltonian. We have then a new form for the eigenvalue 
equation from eq(90) 

V-$ = TEp . (91) 



We obtain in the standard way the integral version of eq(90) namely 

1, : dt 



dx®-p= I dS-j , (92) 
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where R is an arbitrary region of space and S is its bounding surface. 

The norm of the state tp also has a natural interpretation. If we define U so that 

U=\{^^) , (93) 

then we find 

U = \j dxp{x)V{x) . (94) 
That is, the norm of has can be interpreted as the electrostatic energy of the state. 

6.2 One Fermion Level 

A parallel analysis can be carried out at fermionic level one. A state ip e has the form 

%l> = H n (x)a\sl> , (95) 

and is annihilated by Q . That is 

*QV - (d m - lp<j>m(x))ipo = [(d m - P<j> m (x))H m (x)]i> = . (96) 

It follows that 

d m Po(x)H m (x) = . (97) 
That is, the state is associated with a conserved vector field, B n (x), where 

B n (x) = P (x)H n (x) . (98) 

If we adopt eq(44) as the equation of motion, then we have 

® t P (x)H n =Td m [P (x)(d m H n -d n H m )] . (99) 

We can interpret this equation in terms of an electromagnetic analogy in the following way. 
We identify H = (H\, H 2 , H 3 ) as the magnetic field and B = (B\, B 2 , B 3 ) as the magnetic 
induction field. Eq(98) then implies that B = Po^H and hence that Po(x) plays the role of 
(position dependent) magnetic permeability. We introduce an electric current j and an electric 
field E. We link them through a space dependent version of Ohm's Law, 

j = a(.T)E , (100) 

where the space dependent conductivity a(x) is given by 

- ik ■ (101) 

Eq(97) is easily restated as 

V • B = . (102) 
Eq(99) is reproduced by combining Faraday's law of induction 

^B = -VxE , (103) 

with Ohm's law and Ampere's law 

V x H = j . (104) 

That is we have 

-B = -TV x (Pq(x)V x H) . (105) 
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By setting H = H^e TEt etc. where E is an eigenvalue of the supersymmetric Hamiltonian 
at the one fermion level, we obtain the eigenvalue equation in the form 

V x (Po(ic)V x H<°>) = EP (x)U^ . (106) 

It is interesting that allowing for the position dependent magnetic susceptibility represented 
by Pq(x) this equation has a form analogous to the London Equation for a superconducting 
medium [21, 22]. 

We can obtain an integral version of eq(103) in the form 



dt 



[ B • dS = - [ E • dx , (107) 
Js Jc 

where S is an arbitrary surface and C is its bounding loop with the conventional right-handed 
orientation. We also have the integral version of eq(104) namely 

f j • dS = [ H • dx . (108) 
Js Jc 

Again, from the norm of the state, we can identify an energy U given by 

U=\{^) = \ JdxU-B . (109) 
In the context of the electromagnetic analogy this a very natural result. 



6.3 Two and Three Fermion Level 

For completeness we discuss briefly the two and three fermion sectors. 

The evolution of the states in the subspace can be obtained by the application of 

to the states in that we have just discussed. 

The evolution of the states in can most easily be investigated by applying Q to the 

corresponding states in H.^ \ These states, because we are in three dimensions, have the form 

V> = M(x)^ei mn a\al n al l ip . (110) 
The supersymmetric equation of motion then implies 

^M = Tdi{d l -P<f>i)M ■ (HI) 

This is very close to the original diffusion equation for P(x, t) but with the replacement (f> — >• — 4> 
[13]. However since cxp{/3(j)(x)} is not a normalisable wavefunction, there is no stationary 
solution corresponding to a zero eigenvalue of the Hamiltonian H . 

7 Principle of Minimum Dissipation 

We derive here the Principle of Minimum Dissipation in the context of the electromagnetic anal- 
ogy. We will be able to make use of the principle to derive the spatial structure of wavefunctions 
with low lying eigenvalues in the large (3 limit. 

7.1 Zero Fermion Sector 

Eq(90) implies that the energy dissipation rate W is 



W 



= = J dxVd n j n = - J dx(d n Vj n ) . (112) 
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The last step is achieved by integrating by parts. The result is clearly consistent with the idea 
of an electric current being driven by by a potential gradient. It follows that 

w=p l dx A) ■ (m) 

It is a well known principle that, for a given set of sources, the true current distribution 
minimises the energy dissipation as calculated from eq(113). We can verify this in the present 
context by noting that if we modify j n so that 

j n = -TP (x)d n V + d m A mn , (114) 

where A mn is an arbitrary antisymmetric tensor field then d n j n is unchanged and therefore 
j n retains the same (time-dependent) source distribution as the original current distribution. 
However the new dissipation rate calculated from eq(113) is 

W = J dxTP (x)(d n V) 2 + J dxj^(d m A mn ) 2 . (115) 

The cross term vanishes after integration by parts. Clearly the minimum of W occurs when 
d m A mn vanishes and j n acquires its physical value. This outcome can be described by saying 
that for given sources the true current follows the path of least resistance. Under appropriate 
circumstances we can use the principle to identify (or guess) the true current distribution given 
the source distribution. For example we will be interested later in circumstances in which is 
very large. In this limiting regime the hills and valleys of 4>{x) are highly magnified in the rise 
and fall of the local resistivity, (Pq(x))^ 1 . Indeed certain regions where (f)(x) is large become 
effectively insulators. We are therefore in a good position to identify the complementary regions 
in which the current can flow easily and in which in which j 2 remains substantial. The structure 
of these regions depends on the detailed shape of 4>(x) but in the models we consider they will 
be filaments along lines of steepest descent. The calculation becomes essentially one in electric 
circuit theory where the ideas of minimal dissipation are very familiar. 

Of course approximation schemes based on the principle of minimal dissipation could be for- 
mulated to deal with circumstances in which the temperature is well above the low temperature 
limit. Wc do not pursue such calculations in this paper. 



7.2 One Fermion Sector 

Equivalent results can be obtained at the one fermion level. From eq(105) we can verify that 
the dissipation W is given by 



W = = J dxP Q {x){V x H) 2 = (3 J dxj^E 2 = J dxj ■ E . 



(116) 



We can modify E by making the replacement E — > E new = E + WA where A is a scalar field. 
This leaves the law of induction unchanged. The dissipation however is changed to 

w = » I d ^ E "" )2 = !i S ix T^)*' W^™'- (117) 

It is clear that the minimum occurs when \7A vanishes and E acquires its true value. We can 
see from this principle that for a given set of sources the electric field will be greatest where 
Po(x) is greatest and least where it is least. In the low temperature limit these variations are 
exaggerated and this allows us to separate the space into regions of essentially zero electric field 
and complementary regions, confined to thin sheets in our model, where E is large. 

Of course these results can be generalised to any fermionic level but we shall not present 
them explicitly. They will help us in understanding the structure of the eigenstates of the 
supersymmetric Hamiltonian and their associated current distributions. 
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8 High-,5 Limit and Low Lying Eigenvalues 



Witten [14] emphasised the importance of the limit of high j3 in analysing the structure of sta- 
tionary points of 4>{x) . Tanase-Nicola and Kurchan [15, 16] interpreted the limit in the context 
of diffusion processes and developed computational algorithms for exploring the neighbourhoods 
of stationary points. The model they considered was one in which 4>{x) is a Morse function and 
therefore has stationary points with non-degenerate quadratic structure. We discuss the model 
here and show how we can develop the electromagnetic analogy to compute decay exponents in 
the high /3-limit. 



8.1 Zero Fermion Level 

We will investigate a model in which 4>{x) is a Morse function with a set of R minima at the 
(sufficiently well separated) points x — x^ where r = 1, . . . , R . The minima have a quadratic 
structure with the result that in the neighbourhood of the minimum at x^ we have 

4>{x) =B r + l -{x m - x$)A$ n {x n - i«) , (118) 
where B r = (j){x^) and Amn = <j> mn {^ r ^ ■ Clearly 

A(r) _ \( r ) c ( r ) c ( r ) 
7 

where are the eigenvalues of </w(5^ r ') and c^n are the corresponding eigenvectors. We 
are interested in a situation in which (3 is very large (low temperature). In these circumstances 
cxp{— (3(f)(x)} is significantly non-zero only in the neighbourhoods of the minima. We have 

exp{-^(x)}^2exp{-/3[B r + l(x m -xW))^),( a;n -4'-))]} . (120) 

r 

This approximation is adequate for computing the limiting expectation values of appropriately 
smooth and slowly varying functions of x . However we note parenthetically that there are 
circumstances in which we may wish to compute the next-to-leading corrections and that this can 
be achieved by including terms in the eigenfunction approximation that involve third derivatives 
of 4>(x) . The effect can be accounted for by shifting the mean of x m from by an amount 
that is 0{I3- V ) . 

The partition function, in this approximation, can be expressed in the form 

Z = Y,^' 3Br Z {r) , (121) 

r 

where 

Z (r) = J dxc W {-^m - x^)A^ n (x n - iW)} . (122) 
If we define w r so that 

w r = y/e-PBrZW/Z , (123) 

then we can interpret u>J? as the probability of the particle being in the neighbourhood of the 
rth minimum. The stationary probability distribution is approximately 

Po(s)*5>? P o (r) (*) , (124) 

r 

where Pq^x) is the probability distribution appropriate to the rth minimum. Similarly the 
null state of H can be expressed as 

Mx)-Y, W rPo\ X )\°) ' ( 125 ) 
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where (x) = \J P^ r \x) . 

In fact we can use the states i/> (r) = Pq\x)\Q) as an essentially orthonormal basis for the low 
eigenvalue subspace of Hq . The eigenstates of H in this subspace have the form 

iP a =v a (x) Po (x)\0)~J2 v ^ ir ^ w ^ {r) > ( 126 ) 

r 

where v a (x) is slowly varying as a function of x and vq{x) = 1 . Because of the effective 
orthonormality of the eigenfunctions and of the basis set the matrix with elements v a (x^)w r 
is orthogonal, that is 

^ v a (x^)w r v a >(x^)w r = S aa > , (127) 

r 

and 

^ v a (x^)w r v a (x ( - r ^)w r > — S rr r , (128) 

a 

We can extend the basis for the subspace Hq in two ways. First we can add excited states as- 
sociated with the quadratic structure at each minimum. These states are, like the corresponding 
ground states tjj , peaked about each minimum. They are of course orthogonal to these ground 
states. The associated eigenvalues are proportional to f3 . Second we note that, in the limit of 
large /3, all stationary points, including saddles and maxima, yield minima of the potential in 
the Hamiltonian H . We can include then the ground states and higher states associated with 
the quadratic structure of each of the stationary points. These states, in the high-/? limit, are 
all orthogonal to one another and to the states associated with the minima of <fi(x) discussed 
above. Again all of the states have eigenvalues that are proportional to /3 . Of course we must 
cut off the eigenvalue expansions before the corresponding eigenvalues rise to the point where 
we expect serious mixing between the true eigenstates of H . However if the quantities we wish 
to evaluate are dominated by the lower lying of these large eigenvalue states this mixing will 
not be significant. 

We can account for the lowest of these eigenstates that are peaked at saddles rather than 
minima by extending the range of r to label all stationary points and defining the corresponding 

(r) 

states ipQ ' to be 

4 r) =p£\x)\0) , (129) 

where 

P { o\x) = ^ cxp |-|[B r + \(x m - x$)A$ n (x n - 4 r) )]} , (130) 

and 

Amn — 1^7 ^l C 7ro C 7ra ' (131) 
7 

\ip being the eigenvalues (not all positive) of 4> m n(x^) and c^n the corresponding eigenvectors. 
This form for these eigenstates implies that it is not very convenient to express them in the 
"typical" form of eq(86). Of course the higher excitations based on the quadratic structure 
associated with Amli in eq(131) are also available as eigenstates. 

8.2 Fermionic Levels 

An analysis of low lying states at the first fermionic level can be carried out along the same lines. 
As pointed out by Witten [14] and further discussed by Tanase-Nicola and Kurchan [15], the 
relevant wavefunctions are those concentrated at the stationary points x^ of 4>{x) for which 
4>mn{x^) has one negative eigenvalue, that is at saddle points with a Morse index of unity. 
Indeed they contain as factors the saddle wavefunctions discussed in subsection 8.1. The term 
4>mn{x^)a\ n a n in the supersymmetric Hamiltonian appropriate to this saddle, guarantees that 
the lowest lying state has zero eigenvalue in this approximation. Wavefunctions associated with 
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excitations at these saddles have eigenvalues proportional to f3, The other wavefunctions based 
on the original minima and on saddles with a Morse index of two or higher are also available as 
cigcnfunctions but have lowest eigenvalues proportional to /3 . We will not repeat the analysis 
in detail but mention two points. 

First, since the R—l low lying eigenstates in can be raised by the application of to 
provide the low lying eigenstates in there must be at least R—l saddle points of (j>{x) with 
unit Morse index. If there are more than R—l such saddles then the excess is accounted for 
by the corresponding number of low lying eigenstates in By again applying to these 

states we create a number of low lying eigenstates in %^ which then in turn provide a lower 
bound on the number of saddles of <p(x) with a Morse index of 2. The argument can be carried 
on until we run out of higher saddles or maxima. This argument is the basis on which one can 
derive the Morse inequalities from supersymmetry [14, 15]. 

Second the action of on the approximate form for the low lying eigenstates of, say, 
Hq + ' produces zero because the eigenvalue of each of these states is zero in this approximation. 
It follows that the approximation is inadequate for dealing with the subtle structure of the 
eigenstates that allows to move the peaks in the wavefunctions in at the minima of 

<f>(x) to the peaks of the wavefunctions in at the saddles of 4>{x) . 



8.3 Current Flow at Large Times 

As emphasised by Tanase-Nicola and Kurchan [15, 16] the way to understand the link between 
fermionic levels is to study the currents associated with the states in "Hq + ^ . Any initial probabil- 
ity distribution can be expressed as a superposition of eigenfunctions of Hq (see eq(12)). After 
a sufficient length of time the influence of contributions from eigenfunctions associated with 
high lying eigenvalues will have died away and the system will acquire a probability distribution 
that lies in the subspace spanned by the set (x) . As time progresses the occupation of the 
various minima will settle down to the equilibrium values. The equilibriating process is the flow 
between minima of the currents mentioned above [15, 16]. For our purposes then the system 
starts in a state tp where 

ij = V(x)p (x)~J2VrW r p { r) (x) , (132) 

r 

where V r = V(x^) . This corresponds to the probability distribution 

P(x) = V(x)P (x)~J2 V r W r P r \x) ■ (133) 

r 

We assume that V(x) is, in the neighbourhoods of each of the minima, a relatively smooth 
function after the initial fast partial relaxation. Now by exploiting our electromagnetic analogy 
we can identify a charge, q r , associated with the minimum at x^ r \ 

q r = J dxV r w 2 r p( r) (x) =w 2 r V r , (134) 

and we can identify V r as the voltage at x^ r \ then we see that w 2 is the capacitance at that 
minimum. Now 

tI = -E j «. ( 135 ) 

s 

where I rs is the total current flowing from x^ to x^ . Because of the linearity of the diffusion 
process we have 

I rs =R; s \V r -V s ) , (136) 

where R rs is the resistance between the two minima. As a result the time dependence of the 
equations for the voltages becomes 

™f^r = -HRrs\V r -V s ) . (137) 
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We can then find the low lying eigenvalues of H by looking for the solutions of this equation 
that have the form V r = U r e~ pt where U r are constants that obey the homogeneous equations 

W 2 rP U r = Y,Rrs{Ur-U s ) . (138) 

Of course p = TE, where the values of E that make these equations soluble, are the low-lying 
eigenvalues of H . 

The problem then reduces to calculating the resistances. Recall that the principle of min- 
imum energy dissipation tells us that the actual current distribution is along paths of least 
resistance. The local resistivity is (3/Pq{x) and is least in the neighbourhoods of paths of steep- 
est descent (where these exist) passing between minima and over saddles (Morse index 1). When 
f3 is very large any departure from these neighbourhoods penalises the flow of current in terms 
of energy dissipation very severely. 

Consider a typical pair of minima connected across a saddle at x$ ■ Eq(90) tells us that 
away from the sources of current at the minima the current distribution is divergenceless. This 
confirms that the total current flowing parallel to the curve of steepest descent is constant along 
the curve. Near a point xl on the curve of steepest descent with unit tangent vector t the 
current distribution is j(x) oc t . By definition 

j = -TP (x)VV(x) . (139) 

This tells us that VV(x) oc t near the curve of steepest descent and hence that V(x) is constant 
across the profile of the current. The total current is / where 



1 = J dxu-t , (140) 

x±_ are a set of two Cartesian coordinates orthogonal to the curve at xl and dx± is the transverse 
area clement. 
Now 

P (x) = Pt>(x L )exp{~PxlA ± (x L )x±} , (141) 

for points near x^ . Here Aj_(xi) is the 2x2 matrix determining the locally quadratic cross- 
section of the "wire" conducting the current in the neighbourhood of the curve of steepest 
descent from the saddle. We can define a cross-section value for the wire, A(xl) through the 
equation 

J dx ± P (x) ~ P {x L )A{x L ) , (142) 

That is 

(I) - (j ) ^ • (M3) 

where 71 and 72 are the positive eigenvalues of the matrix A(xl) = {4>mn{xL)} and hence also 
of Aj_(x L ) . It follows then that 

I = -TP Q (x L )A(x L )t -VV(x L ) , (144) 

or 

t-vr(«i)-|vM--/. ft( J^y ■ ("») 

where s is distance along the curve of steepest descent. On integrating along the curve of 
steepest descent between the minima we see that the resistance to the flow of current is 



R = P ds . ... , , (146) 
P (x L )A(x L ) 



18 



where the integration covers the path between the minima. However in the limit of large f3 
the integral is dominated by the contribution from near the saddle at the point x$ where the 
integrand is largest. We have then in this limit 



R = & f dse- 1 ^ 2 = ^ J— (147) 

Po{xs)A(xs)J Po(x s )A(x s )p lll ' [ ^'> 

where —711 is the negative eigenvalue of the matrix A(xg) with components 4> m n{xs) ■ In the 
same way we can calculate the resistance between any two minima of <p{x) . 

Later we will require the above result specialised to one and two dimensions. The one- 
dimensional result is obtained by setting A(xl) = 1 and we have 



"-P7-J-- < 148) 

Po(xs) V 7|| 

The two-dimensional result is obtained by setting the one-dimensional cross-section A(xl) to 



A(.r L ) = x l^ , (149) 



where 71 is the single positive eigenvalue of A(xs) ■ We have 



R=jJ!—,\f^ (150) 
Po(x,s) V 7|| 



8.4 One-Fermion Level 

The behaviour of the system at the one-fermion level parallels that at the zero-fermion level. As 
discussed in subsection(6.2) the wave functions in are characterised by a field distribution 
H(x) and an associated field B(x) = P (x)H(x) satisfying V.B(x) = . That part of an initial 
field distribution represented by eigenfunctions with high lying eigenvalues will decay rapidly 
as will the corresponding contributions to the energy U = 1/2 J dxB 2 / Po(x). The remaining 
contributions from the eigenstates with low lying eigenvalues will survive for much longer and 
will be associated with field distributions B that permit the energy U to be as low as possible in 
the circumstances. In the zero-fermion sector this argument implied that the density retreated 
to the neighbourhoods of the minima of <p(x) . In the one-fermion sector, because the field B has 
zero divergence, it cannot be concentrated only in the disjoint neighbourhoods of the minima 
of 4>{x) but must instead be distributed along the filaments surrounding the paths of steepest 
descent joining these minima across saddles with Morse index one. In other words B lies in the 
same channels as the currents we discussed in the zero-fermion sector. Since B is a conserved 
field the actual distribution must be a superposition of loops of field. This is the same picture 
as elucidated in [15, 16] except that in our electromagnetic analogy we replace the conserved 
current with the magnetic induction field. The number of independent loops is of course equal 
to the number of saddles with Morse index two. 

For simplicity consider a situation with only one saddle of Morse index two and one associated 
loop L. We define a magnetic loop integral, M, by 

M = J dx -H , (151) 

where the integration is round the loop L. As argued above, in the limit of large /3, and for 
times sufficiently large, the B field distribution will closely follow the loop of curves of steepest 
descent as will therefore the field H in the neighbourhood of the loop. If we integrate, as in the 
resistance calculation in the above discussion of the zero-fermion sector, over the cross-section 
of the loop current then we find that the total magnetic flux, B, which is constant round the 
loop, is related to the magnetic field at the loop by 

Bt = P (x L )A(x L )H(x L ) , (152) 
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where t is the unit tangent vector at xl on the loop. That is, along the loop, 

H(x L ) = — -, — ^— — rBt , (153) 
V L! P (x L )A(x L ) V ' 

and therefore on integrating round the loop we have 

M = TR L B , (154) 

where is identical to the total resistance round the loop evaluated as the sum of the resistances 
of the various segments evaluated previously. 

To complete the dynamics we apply eq(107) by choosing the surface S to be a disc-like 
element crossing the channel of steepest descent in which the magnetic flux field is concentrated. 
This implies that 

dB 



dt 



E-rfx . (155) 



c 



Note that this is true for any loop C encircling the the steepest descent curve. The result is 
therefore independent of the position x$d at which the loop C cuts the surface of steepest 
descent emerging from the saddle with Morse index 2. It is also independent of the angle at 
which it crosses the surface. That this is possible is due to the fact that E = TP (x)j is 
significant only in the neighbourhood of the surface of steepest descent, j being smooth in this 
neighbourhood. Independence of the angle of incidence requires that E and j = V x H are 
orthogonal to the surface of steepest descent. On evaluating the right side of eq(155), in the 
limit of high (3, we find 

— = -TP (x SD )T(x SD )n.V x H(x SD ) , (156) 
where T{x SD ) is the thickness of the steepest descent surface defined so that 

n.d X P Q (x) = P (x SD )T(x SD ) , (157) 



where n is the normal to the surface of steepest descent at xsd- If we make explicit the quadratic 
structure of the transverse dependence of 4>{x) then 

P (x) =P (x SD )e-^ 3a ^^ 2 , (158) 

and 

T(x SD ) = , (159) 

where a± is the second derivative of 4>{x) in the direction of n. 
On the surface of steepest descent we have from eq(155) 

i7, " w| - iaki i • (160) 

Integrating both sides of this equation over the surface of steepest descent we obtain 

M = -f dS— ~. (161) 

J SD Po(xsd)T(xsd) dt 

In the limit of large (3 the integral on the right is dominated by xs 2 , the saddle with Morse 
index 2. We have 

HhT HhT P dB 

Here a\ and «2 are the eigenvalues of the matrix of second derivatives of 4>{x) along the surface 
at the saddle point. 
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Combining eq(162) and eq(154) we find 
dB 1 



^TR L P {xs 2 )T(x S2 )V^B . (163) 



From this equation we can read off the decay exponent v. It is 



v = ^TR L P (x S2 )T(x S2 y^ = TR lP^s 2 )^^ • (164) 

This result can be generalised to deal with situations involving more than one saddle with Morse 
index two and hence more than one loop. 

The result can be specialised to the two-dimensional case by setting the thickness T(xsd) 
to unity with the result 

v= ^-TR L P {xs 2 )V^^ ■ (165) 

9 External Electric Field 

We can apply our understanding of the supersymmetry structure to the case of an external 
electric field and the calculation of the associated susceptibilities. In particular we derive the 
appropriate form of the Einstein relation between the susceptibility matrix and the fluctuation 
correlation matrix. 

From subsection 3.2 we see that we are dealing with a field of the form f na {X) = d n f a (X) . 
It follows that 

Q f /maWaU=0 . (166) 
From eq(59) and the fact that commutes with the Hamiltonian H we find 

HQ^ a =0 . (167) 

Now Q^ip a € %2 and H has no zero eigenvalues in this subspace it follows that 

Q^ a = , (168) 

and therefore that ip a eR[^ . 
We have also 

fma{x)a\ n ip Q =iQ^ fa{x)ipo ■ (169) 
Combining this result with eq(59) and eq(168) we find that 

Q\Qlpa-^Pfa(x)lb o )=0 . (170) 

It follows that 

Qip a - i/3f a (x)ipo = -iftcatpo , (171) 
for some coefficients {c a } . Taking a scalar product with i[) we obtain the result 

Ca = (V0,/a(zM)) - J dxP (x)f a (x) = f a . (172) 

We have then 

Qlpa=0(fa(x)-f a )% . (173) 

Multiplying by (fb{x) — /&) and taking the scalar product with tp we find 

((/&(*)- fbWo,Qi>a) = m ab , (174) 

where 



-'ab 



((fa(x) ~ fa)(fb(x) ~ f b )) • (175) 
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The left side of eq(174) can be put in the form 

(QH.h(x) - fb)ipo,ipa) = {-id n fa{x)a) n il} () ,F ma (x)al n ip ) . (176) 
The right side of this equation can be expressed as 

i J dxP (x)(d m fa(x))F ma (x) = i(£ ma d m f b (x)) . (177) 

Combining eq(174) and eq(177) we obtain 

(Smad m f b (x)) = 0Hab ■ (178) 

This is our general form for the Einstein relation. In this form it can be generalised to arbitrary, 
possibly compact, manifolds. 

If we specialise to the case for which f a (x) — x a we obtain the more standard result, relevant 
to a uniform external field, 

tab = P((x a -x a ){x b -x b )) . (179) 

The argument given here is a generalisation of the corresponding argument for a one-dimensional 
system presented in [5]. 

9.1 Computation of Susceptibility 

We are particularly interested in the susceptibility of the system in a uniform external field. 
Note that there is no distinction here between the labels {a, b, ...} and {m,n, ...} . To compute 
this susceptibility we evaluate the conditional susceptibility F mn {x) in terms of the eigenstates 
of H that lie in H 



(+) 
o 

Eq(168) implies that 



^ n =iQ X Xn , (180) 

for some \ n e Hq + ^ . If we set \n = V n (x)ipo then 

F mn {x) = d m V n ( X ) . (181) 

Although F mn (x) is unchanged by the addition of a constant to V n (x), in fact V n (x) and hence 
Xn are rendered unique by the requirement that \n is orthogonal to tpo . Eq(173) implies that 

QQhn = P(x n ~ S„)Vo • (182) 

That is 

H X n = P(x n - x n )^ (183) 

Again denote the eigenstates of H in with eigenvalues E a , by v a (x)ipo. We have then 

Xn =^2C an v a (x)ip . (184) 



where 

From eq(183) we have 



C an = (v a 1p 0l Xn) ■ (185) 



C an = ^r(v a ipo,P(x n -x n )ipo) = ^rd an , (186) 
where d an is the dipole moment associated with the state v a (x)ip and is given by 

d an = J dxP (x)v a (x)(x n - x n ) ■ (187) 
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We have then 

V n {x) = '£v a (x)-g-d cm . (188) 
From V n (x) we can compute F mn (x) . We have 

P (x)F mn (x) = V) j£>(z)-|-<U, (189) 

where jm ^ {x) is a current associated with v a (x) and is given by 

j£ } = Po(x)d m v a (x) . (190) 
The eigenvalue equation for v a (x) is equivalent to the result 

d m j£\x) = -E a v a (x)P (x) . (191) 
This can be re-expressed in the form 

j£\x) = E a (x m - x m )v a {x)P a (x) + d k [(x m - x m )j { k a \x)} , (192) 
where the introduction of the term with x m is arbitrary but convenient. We have then 

Po(a;)F mn (a;) = p^(x m - x m )v a (x)P (x)d an + d k Y kmn . (193) 

a 

where 

Ykmn ^ ] ~J^(x m — X m )j k {x)d an . (194) 
a 

On integrating both sides of eq(193) the total divergence drops out and we obtain 

£mn — ^ ^ d arn d an . (195) 
a 

The result is, as it should be, symmetrical in m and n even though this symmetry is not 
required for F mn (x) itself. Note that the summation over a can be extended to include the null 
contribution from the state Vo • 

9.2 Derivation of Susceptibility at High f3 

In eq(195) the summation ^2 a can be split into two parts. The first part which we denote also 
by ^2 a , has a restricted to the low lying eigenstates of H including the ground state -00 ■ If 
we denote this contribution by ^mn then eq(195) tells us that 

a.rr' 

If we exploit the orthogonality of the matrix with elements {w^v a (x^)} we have 

J2^>a(x^)w 2 r ,v a (x^)=w 2 r 6 rr/ , (197) 

a 

and ultimately 

a^E^^-^^'-^) • ( 198 ) 

r 

The second part of the summation is over high lying eigenstates of H n . Wc can describe 
this summation by making the replacement a — > (r, 7) where r runs over the set of minima and 
7 runs over the high lying states of the quadratic approximation to the r-th minimum. In order 
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to compute the local dipole moment at the r* n minimum we can restrict the summation over 
7 taking into account only those states account for the lowest set of such states for which the 
corresponding v a (x) is linear in x . It is then reasonable to assume that the mixing between 
the neighbourhoods of different minima experienced by yet higher lying states can be ignored. 
There is a subtlety relating to the normalisation of the states. We have assumed that 



dxP (x)v a (x)v a r(x) — 6 aa > . (199) 
For states concentrated in the neighbourhood of the r-th minimum this means that 

^xw^PQ T \x)v a (x)v a '(x) — 5 aa > . (200) 

If instead we denote these localised states by (x) with the normalisation 

dxP^ r) (x)v^(x)v { ;?(x)=6 11 > , (201) 



then we must reconcile the two normalisations by making the replacement v a (x) — > w r 1 v^\x) 
in eq(187) . Denoting then the second contribution to the right side of eq(195) we have 

c /3]T w 2 r f dxP^(x)v^(x)(x m - x m ) f dx'P^\x')v^{x'){x' n - x n ) . (202) 

— (r) 

Finally we can make the replacement x n — > x n in this equation without changing the result. 
It follows then that 

4 2 > ^/3]>> r 2 lt> . (203) 

r 

The result is 

Un = P^W?{(x£-X m )(4 ) -Xn) + X$ n } , (204) 
r 

and is consistent with the Einstein relation if we identify the variance matrix as H mn where 

r 

It should be noted however that the terms sti are 0(1//?). In order for the asymptotic 
expansion to be consistent wc should include all effects of this order. These are accounted for 
by reinterpreting it' as the mean position of a particle that is in the neighbourhood of the 
r* n minimum rather than the position of the minimum itself. The two positions differ by terms 
that arc 0(1/(3) . This difference comes about from an approximation to the ground state 
wave function that goes beyond the quadratic approximation to <f>(x) and includes third order 
derivatives of (j>(x) in its Taylor expansion. We will encounter examples of this in section 10 and 
section 11 . With this adjustment the result is indeed intuitively plausible and in this form it 
can be derived directly from the partition function in the high-/? limit. 



9.3 Area Susceptibility 

Evaluating a higher order susceptibility such as the area susceptibility is a more complicated 
task which we do not pursue here in detail. However, it is of interest to investigate the possibility 
of a generalisation of the Einstein relation that holds for the linear susceptibility. We follow a 
line of reasoning for the area susceptibility that is as close as possible to that for the original 
Einstein relation. The result however is not as useful for reasons we explain. The starting point 
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for the analysis is eq(67) . We continue to work with the electric field case so we have from the 
right side of eq(67), 

fma(x)i>b ~ fmb(x)i>a = ((d m fa(x))lpb ~ (d m fb(x))tp a ) (206) 
= iQH(fa(x) - fa^b - (fb(x) - fb)lPa) , 

where we have used the fact that Q^ a = 0. It then follows that 

QHfma(x)lpb- fmb(x)lp a )= . (207) 

Since commutes with the Hamiltonian H we see from eq(67) that 

HQ^ ab = , (208) 

and hence that 

OVa6 = • (209) 

Eq(67) can now be expressed in the form 

QHQ^ab - iP((fa(x) ~ fa)A ~ (fb(x) ~ fb)4>a)) = . (210) 

That is 

Q4>ab-ll3({fa{x)-fa)4>b-{fb(x)-fb)i>a)=Q^ a b , (211) 

where n ab = g a b(x)->Po e . We have 

Hw ab = QQ^ab = -H3Q((fa(x) - f a )^b ~ (fb(x) - fb)%) ■ (212) 

Eq(212) is sufficient to determine ir a b in terms of f a (x) and ip a ■ We take the scalar product of 
each side of eq(211) with 

(fc(x) - fc)QH.f d (x) - j d )ih ■ 

It is simple to see that the left side can be put in the form 

(Q\fc(x) - f c )Q^ (fd{x) - fd)lp0,1pab) = -({£matinb-£mbtina)fm,c{x)fnd(x)) = ~Xab,cd , (213) 

where Xab.cd is the area susceptibility tensor. However when we complete the manipulation 
by forming the appropriate scalar product with the right side of eq(67) we obtain after some 
tedious algebra two contributions. The first is 

\fl{{fkc{x){f d {x) - f d ) - f kd (x)(fc(x) - f c )m b (fa(x) - fa) - tka(.fb(x) - fb))) , 

and the second is 

-\{{dk9ab{x)){{f kc {x)){f a {x) - f a ) - (f kd (x))(f c (x) - f c ))) . 

These expressions are not particularly helpful. We have presented them however in order to 
underline the point that in this case there is no expression for the area susceptibility that 
parallels the Einstein relation for the linear susceptibility and involves only the distribution 
Pq(x) . Instead the two terms above require a knowledge of the complete joint equilibrium 
distribution P(x,£) . This reduces the usefulness of the relation and means that the only way 
to approach the numerical evaluation of the area susceptibility is through the full simulation of 
the slave equations. We show some results later for a specific case in section 11 . 
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10 One-Dimensional Examples 



There are also some general issues that can be addressed particularly easily in one dimension. In 
a previous paper [5] we examined the the case of a double minimum profile and emphasised the 
power-law probability distribution acquired by the slave variable £ at large values. The origin 
of this phenomenon lies in the existence of regions on x-space for which 4>"{x) < [5, 8, 11, 12]. 

Such a power-law distribution raises the possibility that even though the mean of £ exists 
its variance may diverge together with moments higher than the second. Since we are inter- 
ested in comparing numerical results with asymptotic calculations at small T this phenomenon, 
which becomes increasingly severe in this limit, affects how we assess the outcome of our simu- 
lation calculation. We therefore develop a criterion for predicting the appearance of a diverging 
variance, in the one dimensional context. 

In the limit of large time the resulting static distribution P(x,£) satisfies 



d_ f d_ 
dx \ dx 



+ P(x, + 0^ (<i>"(x)£ - 1) P(x, = . (214) 



The support for P(x, £) lies in £ > with the boundary condition P(x, 0) = . If we define 

Q(x,p) as 

Q(x,p)= / d^PP(x,0 • (215) 
Jo 

then we know that Q(x,p) is convergent for < p < p c where for p > p c Q(x,p) is divergent. 
Since Q(x, 1) is convergent we know that p c > 1 . 
From eq(214) we find 

" 1 "' p<f>'(x))Q(x,p)-pl3ct>"(x)Q(x,p)+pf3Q(x,p-l)=0 . (216) 



dx \dx 

If we set Q(x,p) = q(x 7 p)p (x) then eq(216) becomes 

H p q{x,p)=p/3q(x,p-l) , (217) 

where 

H p = -^ + ^(<j>'(x)) 2 + (p-lW(x) . (218) 

For 1 < p we can solve for q(x,p) in terms of q{x,p — 1) provided H p does not have a zero 
eigenvalue. This will remain true for p < p c . The critical value p c is determined by the 
condition that H p has a zero eigenvalue when p = p c . There are straightforward numerical 
techniques such as Sturm sequencing for computing the eigenvalues of H p . By varying p through 
an appropriate range we can determine the value of p that yields a zero eigenvalue. If this critical 
value is less than two then we are in a situation where the variance of the distribution for £ is 
divergent. 

We can see directly how a power law distribution for £ induces a critical value for p. Define 
the probability distribution function D(£) so that 



D(0 = J dxP(x,£) . (219) 

Our claim is that for large £ 

D(Z) * ^_ , (220 ) 
for some constant A and exponent a. If we set 

D&p)= [ dt'pDg) , (221) 
Jo 
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Figure 1: The critical exponent a for the ^-distribution defined in eq(220) . 
then, assuming the above asymptotic behaviour for D(£) we see that 

It follows that the critical value at which £p diverges is p c — a . We can use this result to 
compute a and compare with the results of our simulation. 

10.1 Double Minimum Profile 

The Double Minimum Profile occurs in a model for which 

4>(x) = \(x 2 -a 2 f • (223) 

Clearly <fi(x) has minima at the points x = a and x — —a . We label them r = 1 and r = 2 
respectively. A maximum of <j>(x) lies at the origin. We note that 0(±a) = 0, 0(0) = Aa 4 /4 and 
<j>"{±a) = 2\a 2 , <p"(0) = -\a 2 . In fact <j>"{x) < for -a/^f{2,) <x < a/^/{3) . It follows that 
we expect the slave variable to exhibit a power-law probability distribution at large values of £ . 

10.1.1 Critical Moment for Double Minimum 

For this model the Hamiltonian H p in eq(218) becomes 

H p = + - a 2 ) 2 + (p - \)P\{Z* 2 « 2 ) • (224) 
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Figure 2: The susceptibility x f° r the double minimum profile denned in eq(223) . 



In Fig 1 we show results for the parameter choice A = a = 1 . Using a Sturm sequencing 
technique we determine the values of p for which H p acquires a zero eigenvalue for a range of 
temperature T = 1/(3 . For sufficiently large T the critical value p c is greater than 2, thus 
permitting the variance of £ to exist. However as the temperature is lowered there is then a 
critical value T c below which p c is less than 2, and and for those values of T the variance of £ 
diverges. In the present model T c = 0.473 . 



10.1.2 Simulation for Double Minimum 

The stochastic differential equations for our choice of parameters become 

dx 



dt = -x(x 2 1) + W(t) , (225) 



where (W(t)W(t')) = 2TS(t - t') , and 

^ =-(3a; 2 -!)£ + ! . (226) 



dt 

Using the numerical integration techniques outlined in section 4, we can obtain from a sample 
of N ~ 10 6 particles a measurement of D(£,p) as a function of £ . The case p = 1 allows us to 
extract a and \ = £, the susceptibility, by fitting the asymptotic form in eq(222) to the results 
of the simulation over an appropriate range of £ . 

The variance of x can be evaluated by means of a direct numerical evaluation of the integral 

(x 2 ) = i J dxx 2 e'^ . (227) 

Then using the Einstein relation, \ — £ — P( x2 ), we obtain an evaluation of the the susceptibil- 
ity. In Fig 2 we show that our simulation results compare well with this evaluation, even down 
to temperatures well below T c . 

Of course we could extract £ directly from the simulation data. The results coincide with 
those obtained by the more elaborate fitting procedure outlined above for sufficiently high values 



28 




Temperature T 



Figure 3: The Ratio of full susceptibility to asymptotic form at low T for the double minimum 
profile defined in eq(223) . 



of the temperature. However for T < 0.2 a straight average estimate yields results, indicated in 
Fig 2 by open triangles, that are too low. The elucidation of this inconsistency lies in the manner 
in which the power law tail of the ^-distribution is achieved, namely by the appearance of a 
sequence of large spikes in the evolution of the variable £ as a function of time in the simulation 
[5]. These spikes, which are essential for establishing the power law tail of the ^-distribution, arc 
increasingly prominent and increasingly sparse as the temperature, T, is lowered. Achieving an 
accurate representation of the spikes from the numerical integration of the stochastic differential 
equations is increasingly difficult as T is reduced below T c requiring as it does, a longer and 
longer equilibriation time in order to establish the tail of the ^-distribution. The result is that 
the distribution is underestimated at large £ yielding a value for the susceptibility, y = £ that 
is too low. Our results show however that using the fitting procedure described above, it is 
possible to obtain an estimate of D(£,l), over a range appropriately restricted at small and 
large £, that is sufficiently accurate to permit the extraction of a good estimate of y . 
From cq(227) we find immediately the high j3 limit of the susceptibility to be 

This coincides with the specialisation of eq(205) to one dimension. 
For our choice of parameters 

X=^-l • (229) 

In Fig 3 we show the ratio of the numerical evaluation of y to the above asymptotic form as a 
function of temperature for T < 0.2. It correctly approaches unity at T = with zero slope. 
The correction to the asymptotic form is small for T < 0.02 . As indicated above it is difficult 
to perform the stochastic simulation at these low T values. 
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10.1.3 Decay Rate for Double Minimum 

We can calculate the low T estimate for the decay rate of an asymmetrical initial x-distribution 
using the electromagnetic analogy explained above. This requires a calculation of the resistance 
R between the minima. From eq(148) we find for our model 



R=\1^^ = ^Z**\W\ ■ (230) 



Xa 2 P„(0) V Xa 2 
where the partition function Z is given by 



Z= dxe-^W . (231) 



When (3 is very large we can compute Z from the dominant contributions coming from the 
neighbourhoods of the two minima. Both contributions are equal so 



Z 

We have also 

where w\ = w\ = 1/2 and 



= 2 J dxe-P Xa2{x - a)2 = 2^^2 • ( 232 ) 
Po(x)=w 2 p( 1) (x)+w%Pk 2 \x) , (233) 



P (1 ' 2) (z) = \/— e-^ 2 ^) 2 . (234) 



7T 



We note that <p"(0) = -Xa 2 and that 

Po(0) = 4^ Aa4 • (235) 
Applying the electrical analogy we see that the voltages at the two minima satisfy the equations 

1% ' (236) 



and 



where R is given by eq(230) . 



1 dVo 1 

2li = -RM- V ') ' (237 > 



^Vot- (238) 

It follows that the rate of decay of the voltage difference is proportional to e~ vt where 

4 y/2Xa 2 4 

v = - = e -4^« 239 

R IT 

The eigenvalue of the supersymmetric Hamiltonian is of course f3v . 

In Fig 4 we show the results of the simulation (for the case A = a = 1). They compare well 
with the direct evaluation obtained from the Sturm sequencing evaluation of the lowest non-zero 
eigenvalue of H . For this choice of parameters the asymptotic result becomes 

I/= ^ e -V(4T) _ (240) 

For the temperature range T < 0.2 this formula yields results very close to the simulation and 
numerical evaluation shown in Fig 4. 
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Figure 4: A comparison of the decay exponent v for the double minimum profile defined in eq(223) 
calculated theoretically from the Einstein relation and from a simulation numerically integrating 
the stochastic differential equations, eq(225) . 

10.2 Three Minimum Profile 

As a matter of interest we work through the theoretical calculations for a more complicated 
example for which 4>{x) has three minima and two maxima. We choose 

0(x) = ^x 2 (x 2 -a 2 ) 2 . (241) 

The minima are at the points x = 0, ±a, and the maxima are at x = ±a/\/3. We have the results 
(f>(±a/V3) = Aa 4 /27 and <//'(0) = Aa 4 /2, <f>"(±a) = 2Aa 4 together with 0"(±a/\/3) = -2Aa 4 /3. 

If we label the minima at x = ±o by 1 and 2, and the minimum at x = by then we easily 
find in the high j3 limit 



and 



y ■ (243) 



Hence 

Z = Z 1 + Z 2 + Z =A ] j^ , (244) 

and therefore w\ — w\ — 1/4 and Wq — 1/2 . 

The resistance between the origin and each of the other two minima is R where 

R=—e^. (245) 
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If the voltages at the minima are v\, V 2 and Vq then the electrical analogy tells us that 





1 


Ilk ~ 


~R 


1 dV 2 _ 


1 


4 dt 





(Vi-V ) , (246) 

(V 2 -V ) , (247) 
and 

It follows that equations for the eigencombinations are 

^- t (V 1 + V 2 -2V ) = -^(V 1 +V 2 -2V Q ) , (249) 

and 

§ i (V 1 -V 2 ) = ~(V 1 -V 2 ) . (250) 



The two decay exponents are then 



and 



2Aa 4 e\« 4 

^i=y^e-— , (251) 
z, 2 = e - — . (252) 

V07T 

The susceptibility can be computed from eq(205) to be 

11 Two-Dimensional Example 

We can illustrate some of the results of general theory by means of two-dimensional models. 



11.1 Blind Saddle 

In [15] Tanase-Nicola and Kurchan introduce a simple model that they refer to as a blind saddle. 
It is constructed as a two-dimensional tilted Mexican hat for which 

0(x) = ^\{x 2 + y 2 - a 2 ) 2 - Jx . (254) 

In this case, for sufficiently small J, <fi(x) has one minimum near x = a and y = 0, one saddle 
with Morse index one near x — —a and y — together with a maximum, Morse index two, near 
the origin. All lines of steepest descent leaving the saddle arrive at the the single minimum. 
The uniqueness of the minimum means that there is no opportunity to construct states with 
low non-zero eigenvalues for Hq in the high j3 limit. The interest of the model then, is that in 
the limit of large (3 there are no eigenstates in with exponentially small eigenvalues. The 
lowest eigenstates in are all of 0(/3) . However in there is such a low-lying eigenstate 
which we associate with a magnetic flux loop B passing round the circuit created by the lines 
of steepest descent joining the saddle to the minimum. 

To detect this low lying state we must first create it. Subsequently we let it die away and 
measure the exponential decay of an appropriate observable, namely . This yields the value 
of the decay exponent. 



32 



To create the state we integrate the stochastic differential equation, eq(l) together with the 
associated slave equation, eq(49) in which f na (x) is chosen to have a magnetic character. That 
is 

Qfna(x)c&ip = . (255) 

Using the fact that Q commutes with H, we see from eq(59) that the form of ip a (x) in the limit 
of large time, also has the property that 

Qi> a =0 . (256) 

This guarantees that ip a € %^ and therefore corresponds to a state with a loop of flux as 
described in subsection 8.4 . We can solve eq(59) for ip a (x) by introducing an eigenfunction 
basis, {uj a }, for . 

ty a = 2_^u a p • (257) 

If now we continue the simulation by integrating the associated slave equation with f na (x) set 
to zero, the subsequent time development of ip a is given by 

where th time t is measured from the new starting point. Consider the matrix element 

{fmb{x)ahil)o,il>a) = J dxP (x)f mb (x)F ma (x) = (f m b(x)^ ma ) ■ (259) 

It follows that 

(/ b(x)€ ) = ^ c -TE Afrab{x)a\ n ^,u o ){u ai ,f na {x)a\ l 7p ) 

The long time behaviour of this expectation value is governed by the most slowly varying 
exponential. By measuring this expectation value in the simulation we can obtain the decay 
exponent. At this point a limitation of the method is encountered. A careful analysis, not 
presented here, shows that the scalar product (u) a , f na (x)al l ipo) — C aa V E a , where the coefficient 
C'aa depends on the form of f na (x) . The explicit dependence on the eigenvalue E a removes the 
potentially enhancing denominator on the right in eq(260) and for reasonable choices of f na (x) 
the coefficient C aa is sharply reduced at low temperature. This results in the decoupling of the 
contribution from the lowest eigenstate as the temperature is lowered making it harder to pick 
out the relevant term. This barrier to easy simulation at low temperature is a theme running 
through our calculations and prevents a close comparison with the low temperature asymptotic 
result as we see below. 

In practice we need only a single flow field and we drop reference to the label a from now on. 
We choose f x (x) = 4> y {x) and f y (x) = —<fi x (x) . It is easily checked that f n (x)al l ipo € . 
For the purposes of the simulation we set A = a = 1 in eq(254). We have therefore 

(261) 



<t>{x,y) 


= 4 (2; 


+ y 2 - 


- I) 2 - Jx 


4>x{x,y) 


= x(x 2 


+ y 2 - 


1)-J , 


<f> y {x,y) 


= y(x 2 


+ y 2 - 


1) , 


<Pxx(x,y) 


= (3x 2 


+ y 2 - 


1) , 


^yyix.y) 


= (z 2 4 




1) • 


4>xy{x,y) 


= 2xy 







33 



T 


^sim 


"th 


^asym 


0.4 


0.8813(5) 


0.884(2) 


0.8760 


0.3 


0.6449(5) 


0.651(1) 


0.7230 


0.2 


0.3900(6) 


0.392(1) 


0.4925 


0.18 


0.3411(5) 


0.3386(7) 


0.4334 


0.16 


0.283(1) 


0.2835(7) 


0.3693 


0.10 


0.110(4) 


0.1178(5) 


0.1557 


0.09 


0.084(3) 


0.0946(5) 


0.1206 



Table 1: Comparison of values of the decay index v obtained by simulation (f s j m ), lowest eigenvalue 
in T-6 2 ^ (^th) anc ^ ^ e as y m Ptotic estimate (fasym) • 



In order not to shift the stationary points very far from x = 0, ±1 we set J = 0.02 . These 
stationary points which lie on the line y = . On the left we have the saddle at x = —0.989846 
with 4> xx = 7i = 1.939384, <j> yy = ~7|| = -0.020205, the maximum at x = -0.020008 with 
cj) xx = ai = -0.998799, <j> yy = a 2 = -0.999599, and the minimum at x = 1.009854 with 
4> xx = 2.059415, 4> yy = 0.019805 . Of course (f> xy — at all three stationary points. For 
future reference we note that at the maximum we have 0(max) = 0.250200 and at the saddle 
^(saddle) = 0.019899 . Hence the difference in heights between maximum and saddle is d where 

d = 0(max) - ./.(saddle) = 0.23030106 . (262) 

Combining eq(150) and eq(165) we find the asymptotic estimate for the decay index 

1 /olaiTi - a 
^asym-^^-e 

Using the numbers appropriate to the simulation we have 

^sym = ge- d/T , (264) 

where g = 1.558012 . A comparison between the asymptotic formula and the simulation results 
is shown in Table 1. 

Also recorded is our evaluation of the decay index as = TE where E is the lowest 

eigenvalue of H in the sector . Because we are in two dimensions this is the lowest eigen- 

value of the simple Schrodinger Hamiltonian Hq obtained from Hq by making the replacement 
4>{x) — > ~4>{x) . We use the lattice field theory method of reference [23] to evaluate this eigen- 
value. We see that this evaluation compares favourably with the outcome of the simulation, 
particularly at the higher values of temperature. The estimate from the asymptotic formula 
is less good, though not totally out of line with the other results. Our experience with the 
one dimensional problem suggests that we would require values of temperature somewhat lower 
than 0.1 to see agreement with the other methods of evaluation. However as we have indicated 
above these methods are not easily applied for temperatures in this range. With these provisos 
we feel that the theory and simulation are in encouraging agreement. It would be interesting 
to attempt theoretical evaluations of v at the higher temperatures based on an approximation 
scheme that exploits the application of the minimum dissipation principle to appropriately cho- 
sen field distributions that are less extreme in form than the very narrow wires and surfaces 
used in the very low T limit. 



(263) 
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11.2 Four Minimum Model 

Finally we illustrate our methods with a model that has four minima. We choose cj> have the 
form 

= 4>{x, y) = \\(x 2 a 2 ) 2 + - a 2 ) 2 + \nx 2 y 2 , (265) 

where x and y are the coordinates of the diffusing particle. The relevant derivatives are 

4> x — \x(x 2 — a 2 ) + nxy 2 , (266) 

4> v = ny(y 2 - « 2 ) + Kx2 y , 

<\>xx = M^x 2 - a 2 ) + ny 2 , 
4> yy = n(3y 2 - a 2 ) + kx 2 . 
4> xy = 2nxy 



The stationary points satisfy 



x[X(x 2 - a 2 ) + ny 2 } = , (267) 
y[^{y 2 - a 2 ) + kx 2 ] = . 



We have then the following stationary points. 

(1) A maximum at x = y = with <p xx = — Aa 2 , <j> yy = — tta 2 , cf) xy = . 

(2) Saddles with Morse index 1 at x — 0, y — ±o with <f> xx = — (A — n)a 2 , <fi yy = 2[ia 2 , (f> xy — 
and at y = 0, .t = ±a with (/> xx = 2Aa 2 , <f> yy = —(jj, — n)a 2 , <j) xy = . 

(3) Four minima at the points x = ±\J n{\ — k)/(A// — n 2 )a and y = ±^X(fj, — k)/(A/x — n 2 )a 
with <j> xx = 2a 2 Xfi(X — k)/(X^i — k 2 ) and <p yy — 2a 2 X^(n — k)/(A/z — k 2 ) and XJ/ = 
±2a 2 K^/Xn(X — k)(h — n)/(Xfi - k 2 ) . 

For case (3) the discriminant of the Hessian at each minimum is A where 

A = 4a M\- -)(P - -) . (268) 
A/i — K 

We see that these stationary points are indeed minima if < n < min(A, /z) since this guarantees 
that A, cj) xx and <p yy are all positive. The eigenvalues of the Hessian matrix are 

2 

— - — 2 | A M(A + M - 2k) ± VA 2 A* 2 (A - /i) 2 + 4k 2 A/i(A - /t%u - is)} 
11.2.1 Decay Indices 

We first calculate the decay indices in the high fi- limit. In this limit the partition function is 
^ Stt 1 f /3a 4 «(2A//- K (A + /i)) ) 

ex Pi — 5 r- — -5 r • ( 269 ) 



/3 7A I 4 A M 

To calculate the low lying eigenvalues in this limit we first evaluate the resistances between 
adjacent minima. In the present model there are two resistances namely, R x , appropriate to 
paths of steepest descent passing in the ^-direction through the saddles at x — 0, y — ±a, and 
Ry, appropriate to paths passing in the y-direction through the saddles at y = 0, x — ±a . 
Using the result in eq(150) we find 



8tt 



. / 2/i f /3a 4 fi(X - k) 2 . 
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and 

Stt / 9X r Rn* \(n - k\ 2 ~\ 

(271) 

For convenience we label the minima with r = 1, 2, 3, 4 in a clockwise direction starting in the 
positive quadrant. The symmetry of the model is sufficient to ensure that the capacitances at 
all four minima are equal. That is w^. = 1/4 . From eq(137) we obtain the equations of current 
flow as 



IdVi 
4~dt 
l dV 2 
4 dt 

im 

4 dt 

I dVj 

4~dt 



= -(Vi 


- V 2 )/Ry 


-(Vi 


- V 4 )/R x 


= -(v 2 


~ V 3 )/Ry 


-(v 2 


- V 1 )/R x 


= -(v 3 


- V^/Ry 


-(v 3 


- V 2 )/R x 


= -{Vi 


- V 3 )/Ry 


-(Vi- 


- V x ))/R x 



(272) 



The three eigenmodes of the form V r = U r e pt with p^O are 

U 1 = U 2 = U U 3 = U A = -U p = 8/R x , 

U 1 = U i = U U 2 = U 3 = -U p = 8/Ry , (273) 
U 1 = U 3 = U U 2 = U 4 = -U p= (8/R x + 8/Ry) . 

We conclude that the three decay exponents associated with the subspaces and Ji^ ^ 
are v\ = 8/R x , v 2 = 8/R y and v 3 = (8/R x + 8/Ry) . Note that these results at low T imply 
that v 3 = u\ + v 2 . However we do not expect this simple relation to hold for general values of 
the temperature. 

To evaluate these decay exponents for values of T above the low temperature limit we 
performed a simulation by numerically integrating the stochastic equations of motion. Our 
procedure was to cquilibriatc the system and then to fold over the resulting configuration by 
reflecting the position of each particle in appropriate axes so that it ended up in the positive 
(x, y)-quadrant, the remaining three quadrants then being empty. We then restarted the evo- 
lution of the system and measured the expectation values (x), (y) and (xy) as a function of 
time. The long time behaviour of each of these three quantities yields respectively the decay 
exponents u\, v 2 and v 3 . In Table 2 we show the results in the case a — 1.0, A = 1.0 and fi = 0.5 
for a number of values of k and T . The case k = 0.0 corresponds to a situation in which motion 
in the x and y directions are independent. This serves as a check on our procedures since under 
these circumstances we expect that v 3 = v\ + v 2 at all temperatures. Within errors Table2 shows 
that this is indeed the case. For k ^ 0.0 we expect this outcome only at low temperatures, as 
indicated above. At the temperatures listed in Table 2 the result does not hold as indeed we 
anticipated. 

General theory tells us that the dimension of the low lying subspace in "Hi is equal to the 
number of saddles. As we have seen three of the basis states are in so a remaining low lying 
state lies in with a corresponding supersymmetric partner in H 2 ' • In our electromagnetic 
analogy this state is associated with a loop of magnetic flux (in the present example confined 
to the two-dimensional surface) that runs through all the minima and surrounds the maximum 
of (f> at the point (0, 0) . 

The loop resistance is given by 

R L = 2(R X + R y ) (274) 
From eq(165) we can conclude that the associated decay exponent is 




(275) 
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K 


T 


v\ 


V2 


vz 


0.0 


0.8 


0.6573(5) 


0.5704(6) 


1.222(5) 




0.6 


0.5065(4) 


0.4626(6) 


0.9632(2) 


0.2 


0.8 


0.7598(6) 


0.6612(4) 


1.566(3) 




0.6 


0.5847(3) 


0.5452(3) 


1.239(2) 


0.3 


0.8 


0.8000(7) 


0.6985(4) 


1.732(4) 




0.6 


0.6174(3) 


0.5792(4) 


1.373(4) 



Table 2: Simulation results for the decay indices vi, and v% for the Four Minimum model. 



K 


T 




^4(th) 


0.0 


0.8 


1.354(1) 


1.240(2) 




0.6 


0.954(3) 


0.965(2) 


0.2 


0.8 


1.313(2) 


1.317(2)) 




0.6 


1.031(1) 


1.039(2)) 


0.3 


0.8 


1.354(1) 


1.347(2) 




0.6 


1.068(1) 


1.065(1) 



Table 3: Simulation results for the decay index v± for the Four Minimum model together with a 
theoretical evaluation ■ 

The decay index v± is associated with the subspaces n[ + ^ and %^ . We evaluate it at 
higher temperatures by using the same procedure mutatis mutandis as for the blind saddle. 
The results are shown in Table3 . Also shown for comparison is a theoretical evaluation of the 
index, z^th)' obtained using the Schrodinger problem with <p{x) — > — (f>(x) as was done for the 
blind saddle. The two sets of results compare reasonably well providing good evidence that the 
numerical methods work reliably. It is worth noting that for the case k = 0.0 we v z — as 
expected. 

For practical reasons both the simulation and the theoretical Schrodinger calculation of 
^4 cannot be carried out at a temperature sufficiently low that a comparison with the the 
asymptotic calculation is possible. 



11.2.2 Susceptibilities 

We first calculate the susceptibility at low temperature. From eq(205) we see that there are 
two types of contribution from each of the minima. The overall mean position of the particle is 
the origin in this model and the capacity of each minimum is 1/4. We obtain from minimum 1 
contributions to T, xx , S yy and Y< xy of the form ^(x^) 2 , \{y^) 2 and ^(x^)y^ where x^ and 
j/W refer to the position of the minimum modified by terms that are 0(l//3) . This modification 
is computed by setting 

(276) 

, ... - (277) 

\[l — K z 
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We then consider x' and y' to be small quantities and evaluate <j){x,y) in the neighbourhood of 
the minimum as 

4>{x, y) — (^(minimum) + quadratic . (278) 

where 



quadratic = -{<j) xx x 12 + (j) yy y' 2 + 2<\> xy x'y') + cubic , 

cubic = ^{(f) xxx x /3 + <j) yyy y /3 + 3<j) xxy x /2 y / + 3x'y /2 ) . 



(279) 
(280) 



If we retain this cubic term in <p in evaluating the local average of x and y at minimum 1 we 
find to 0(l//3) 



?(!) - 



= (l\ 



'MA 



A/i — K 2 



y 



(i) 



'a(m- 



A/i — K 2 



1 3AV(^-k) + 2k 3 (A-k) 
4,3a 4 A^(A-k) 2 (^-k) 

1 3A/i 2 (A - k) + 2k 3 - k) 
4/3a 4 A//(A- k) 2 



(281) 
(282) 



The average positions at the other minima are obtained by applying the same modifying factors. 

The second set of contributions from eq(205) are obtained by forming the local averages at 
minimum 1 of x' 2 , y' 2 and x'y' keeping only the quadratic approximation for . We find 



„/2 



,/2 



1 



1 



x'y' 



2f3a 2 \-n ' 

1 1 
2(3a 2 fi — k 

1 K 



2(3a 2 y>\^(\- n) 



(283) 
(284) 
(285) 



The the other minima yield identical results except for a sign change in the x'y' from minima 
2 and 3 . The outcome from adding all the contributions from each of the minima is that T, xy 
vanishes and 



yy 



a 



a: 



MA - «) 

A/i — K 2 

2 A(/x-k) 

A/i — K 2 



1 - 



1 



1 2A 2 fi(fi - k) + k 2 (X/i - k 2 + k{\ - k)) 
213a 2 \fi(\- k) 2 (li- k) 

1 2A/i 2 (A — k) + K 2 (A/i — k 2 + k(/j, — k)) 



2(3a 2 



A/i(A — n)(p — k) 2 



(286) 
(287) 



We show in Fig 5 this low-T asymptotic dependence of the susceptibilities for a = 1.0, A = 1.0 
and /i = 0.5 . For the same parameter values we can obtain the susceptibilities from a direct 
simulation of the stochastic differential equations. Guided by our experience with the one- 
dimensional problem we extract the results by fitting an appropriate power law form to the 
measured cumulative distributions for £ xx and £ yy . The results for the susceptibility in the 
x-direction are exhibited in Fig 5 in comparison with those predicted using the Einstein relation 
from an evaluation by direct numerical integration of (a; 2 ) and the asymptotic low T results 
from cq(286). 

There is good agreement between the simulation and the direct numerical valuation for 
temperatures T > 0.15 . For lower values of T here is difficulty in achieving good simulation 
results for the same reason as in the one dimensional case, namely that the distributions for £ xx 
and (; yy acquire such slowly descending power law tails that the do not have finite variances. 
However the asymptotic estimate does agree well with the direct numerical evaluation at low 
values of T . 

In Table 4 we exhibit the results from the simulation for the indices of the power law tails 
at a number of values for k and the temperature T . The first point of note is that for the case 
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Figure 5: The susceptibility Xx = ^xx for the Four Minimum Model calculated in the low temper- 
ature limit, exactly from the the Einstein relation and from the stochastic and slave equations. 

k = and at both temperatures, the area index a xAy is reasonably coincident with the smaller 
of the two linear indices namely a x . The explanation for this lies in the fact that when k = 
the variables £ xx and £, yy evolve independently and the area element is just z = £xx£,yy ■ In the 
appendix we show that a product variable such, as z, has in general, a power law tail with an 
index that coincides with the smaller of the two indices involved in the product. Of course when 
k is no longer zero we cannot infer any such simple connection and the results in Table 4 bear 
this out. 

In Table 4 we have included the area susceptibility for the relevant values of k and T . We 
found that obtaining a good outcome for the simulation at much lower temperatures rather 
difficult given the size, already considerable, of the statistical sample. Although we can in 
principle compute the low T limit of the area susceptibility by generalising previous calculations 
for the linear susceptibility we do not pursue this point here. Since it is not possible to construct 
an equivalent to the direct numerical evaluation the linear susceptibility through the application 
of the Einstein relation we cannot carry out a numerical comparison. Nevertheless an exploration 
of such higher susceptibilities remains an interesting topic of further study. 



12 Conclusions 

We have studied a stochastic system comprising a particle with degrees of freedom {X n } together 
with associated slave variables {S ma } that represent infinitesimal line elements carried along 
by the diffusing system. We refer to the equations satisfied by these slave variables as slave 
equations. The states of the system are described by a joint probability distribution P(x, £), 
which may depend on time and which satisfies an appropriately generalised diffusion equation. 
Here {x n } and {£ ma } are values attained by the stochastic variables {X n } and {S ma } . The slave 
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K 


T 


a x 


a y 




Xarea 


0.0 


0.8 
0.6 


2.76(1) 
2.31(1) 


4.7(3) 
3.8(1) 


2.88(2) 
2.33(2) 


1.835750(3) 
2.8086(1) 


0.2 


0.8 
0.6 


3.62(3) 
2.99(2) 


6.2(2) 
4.79(2) 


3.29(8) 
2.60(1) 


1.5734(2) 
2.38711(6) 


0.3 


0.8 
0.6 


4.07(7) 
3.29(4) 


6.9(6) 
5.64(9) 


3.5(1) 
2.72(2) 


1.47233(2) 
2.22607(4) 



Table 4: Simulation results for the power law tail index for the x-susceptibility (cx x ), ^-susceptibility 
(a y ),the area susceptibility (a X Ay) and the value of the area susceptibility 

variables can be used to construct a hierarchy of infinitesimal areas, volumes etc, that evolve 
either from an initial set of (infinitesimal) displacements of the initial conditions of the stochastic 
system or through alterations of trajectory induced by external forces. In the former situation 
we are concerned with decay exponents that describe the return of the system to equilibrium. 
In the latter we are interested in the (infinitesimal) response of the equilibrium distribution for 
the system to the presence of external forces. This response can be described by a hierarchy of 
susceptibilities, linear, area etc. There are exponents associated with these susceptibilities as a 
result of the fact that the probability distributions of the associated variables have power-law 
tails. These exponents play a significant role in extracting the values of the susceptibilities from 
the simulation data. 

The state of the system can alternatively be described in terms of appropriately antisym- 
metrized moments of the slave variables {£ ma } • In turn these moments may be regarded as the 
wavefunctions in a hierarchy of states that are the basis of a supersymmetric description of the 
original system. The behaviour of the system is controlled by a supersymmetric Hamiltonian as 
explained by [15, 16] . The decay exponents characterising the return of the system to equilib- 
rium are, up to a factor of the temperature T, the eigenvalues of the supersymmetric Hamiltonian 
[15, 16, 17]. The susceptibilities are expectation values of antisymmetrized moments of the slave 
variables and are therefore related in a natural way to the hierarchy of supersymmetric states. 
This supersymmetric structure is the same as that introduced by Witten[14] and investigated 
further by Tanase-Nicola and Kurchan [15, 16] . 

In order to obtain the decay exponents we follow a two pronged approach, theoretical on the 
one hand and computational on the other. In the theoretical approach we exploit an electromag- 
netic analogy and formulate a Principle of Minimum Dissipation that provides a natural basis 
for computing decay exponents. It is particularly effective in the limit of very low temperature. 
It might also be made the basis, in the style of the Rayleigh-Ritz method, for an approximate 
numerical computation of exponents at higher temperatures, though we have not pursued this 
idea here. 

The computational approach we use is to some extent dependent on the number of degrees 
of freedom of the system and for simplicity we restrict our calculations to one and two degrees 
of freedom although there is no barrier to working with higher dimensional systems. For one 
dimensional systems we perform a direct numerical integration of the stochastic equations of 
motion starting with a non-equilibrium distribution of the variable. By measuring the expec- 
tation value of an appropriate observable we can observe its ultimately exponential decay with 
time and measure the decay exponent directly. This works well above a certain temperature. 
Unfortunately at its lower end this temperature range is just marginally above the range of 
applicability of the low temperature theoretical results. However in the one dimensional case 
it is simple to solve the relevant supersymmetric eigenvalue problem directly by the method of 
Sturm sequencing. This produces accurate results over a temperature range that covers both 
the low temperature asymptotic range and the higher simulation range. We obtain good con- 
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sistency in both cases. This confirms the efficacy of the numerical simulation the correctness of 
the low temperature evaluation. 

For systems with two degrees of freedom we again carry out a simulation by integrating 
the stochastic differential equations including the slave equations. Again by measuring the 
expectation values of appropriate variables we can measure the exponential decay that yields 
the low lying eigenvalues of the supersymmetric Hamiltonian. It should be noted however that in 
order to obtain the full complement of these low lying eigenvalues it is necessary to include slave 
variables that respond to behaviour in the higher fermion sector 1-6^ . This is particularly true 
for the case of the blind saddle [15, 16] where the only low lying eigenvalue of the supersymmetric 
Hamiltonian is in this sector. For this model and the four minimum model it is not easy to 
compare directly with the low temperature asymptotic calculations. However qualitatively the 
results are not bad. By exploiting lattice field theory techniques [23] we were able to compute the 
eigenvalue appropriate to the subspace which because of the supersymmetry inherent in 

the problem is the same as the eigenvalue for the Hi subspace. This yielded results consistent 
with the numerical simulation. 

In studying susceptibilities we again use two approaches one theoretical covering low tem- 
peratures, the other numerical applicable at higher temperatures. The theoretical calculation is 
exhibited in subsection 9.2. It shows that at low temperatures the linear susceptibility has two 
contributions, one associated with the positions and occupation of the neighbourhoods of the 
various minima and the other associated with a sum over the intrinsic susceptibilities of each of 
the minima. The numerical approach involves integrating the stochastic equations including the 
slave equations. After equilibriation we can in principle compute the mean value of S ma directly 
from a sample of results obtained by continuing the integration procedure. At sufficiently high 
temperatures this works well. However at lower temperatures, as explained in [5] and subsec- 
tion 10, the probability distribution for the slave variable, £, which has a power-law tail at high 
values of £ . At sufficiently low temperatures the exponent of this power is so small that it 
prevents the slave variable having a variance. This makes the estimation of the average from 
the simulation impossible. However we were able to obtain good results by compiling a running 
average, which does have a variance, and by fitting the appropriate asymptotic form to the 
output from the numerical simulation. We were able, for the one dimensional case, to formulate 
a method for computing the exponent and hence deduce the critical temperature below which 
the variance diverges. This amounts to constructing a modified Hamiltonian H p , associated 
with the moment £ p and then adjusting the value of p until H p acquires a zero eigenvalue. This 
procedure is easily carried out numerically. In two dimensions we obtained similar results but 
without a method of establishing the critical temperature theoretically. This emphasises the 
importance of the simulations approach for more complicated dynamical systems. It is also pos- 
sible to compute the linear susceptibility in any number of dimensions by computing directly 
the variance of the position variables X n by using the Einstein relation (see section 9). We 
obtained good agreement between the different numerical approaches and where relevant with 
the theoretical calculation. In two dimensions it is possible to compute an area susceptibility. 
We were able to obtain numerical results at sufficiently high temperatures from the numerical 
simulation. However we were not able to establish a useful version of the Einstein relation in 
this case in order to provide a check. Further work is also needed to investigate the theoretical 
results for the low temperature limit. 

A question of interest is how the various susceptibilities might be measured. In a simulation 
there is no difficulty since the various relevant quantities can be computed in a direct manner by 
integrating the stochastic differential equation. In fact a measurement of the various correlators 
and decay exponents and their inter-relationships provides a useful means of establishing the 
consistency of the simulation and aids detective work on the structure of the landscape potential 
function. The experimental measurement of the standard susceptibility is, conceptually at least, 
also straightforward. One is looking at the response of the mean displacement to an external 
field. An experimental measurement of higher order susceptibilities remains for the moment 
elusive. However an example, in a different context, in which higher order elements, tetrads 
in fact, play a significant role is in the computation of the magnetic a-effect for magnetic field 
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generation [7]. Finally it is worth pointing out that elucidating the structure and evaluating 
the spectrum of supersymmetric quantum mechanical systems is an intrinsically interesting and 
important goal. 



Appendix 

We analyse here by means of a simple example the asymptotic property of the distribution of a 
variable z that is the product of two independently distributed variables u and v . Because we 
are concerned with distributions on the range < u, v < oo that exhibit power law behaviour 
for large values of the variables we choose as the distributions P(u, a) and P(v, fi) where 

= (° - *) (1 + g 2)(a + l) /a ( 288 ) 

The distribution for z is 

Q(z) = J dudvS(z-uv)P(u,a)P(v,j) . (290) 

The Mcllin transform of P(u, a) is 

/•OO J 

P(l,a)= I -r^P(u,a) , (291) 



and the inverse relation is 



P(u,a) = J dlu l P(l,a) , (292) 

and similarly for the other distributions. The ^-integration is along a contour in the imaginary 
direction with a real part appropriately positioned. By splitting the w-integration range into 
< u < 1 and 1 < u < oo we find explicitly 

oo +1 1 °° +1 1 — 1 

PC,a)-(a- 1) j:c^^ Y - l + { a- 1 )Y.C^^ Ti+ ^ . (293) 

n— n— 1 

The pole at I — —a arises from the behaviour P(u, a) ~ (a — l)/x a as x — > oo . The Z-contour 
passes between this pole and that at I — 1 which controls the behaviour at small u . Similarly 
P(l, (3) has a pole at I = —(3 . 
It is easily seen that 

Q(l) = P(l,a)P(l,i) . (294) 

Therefore Q(l) has poles at I — —a and I — —(3 . If /3 > a then the large z asymptotic behaviour 
will be controlled by the rightmost pole at I = — a . That is 

Q(z)cx^ , (295) 

for large z . The only change to this result arises when /? = a . In this case Q(l) has a double 
pole at I = a and the dominant behaviour at large z becomes 

■ (296) 

Of course if f3 ~ a then it is probably more useful to include both asymptotic contributions. 
However the basic result is that the distribution for z has the same inverse power behaviour as 
that of the most slowly decreasing of the two distributions. 
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